clear all;
clc;
tstart=1; %start time
tmax=50; %max time of simulation
time=tstart:tmax; %time vector

n_iterations=1;
npos=2*length(time)+3; %size of pos arrays = 2t+3
pos=1:npos;
nmid_rel=tmax+2; %relative mid position in the below a arrays which corresponds to zero position in the absolute sense
nstart_rel=0;%relative start position in the below a arrays which corresponds to least -ve position in the absolute sense
nend_rel=0;%relative end position in the below a arrays which corresponds to max +ve position in the absolute sense

a_n_s1_prev=zeros(npos,npos); %a's for spin up at each N at previous time
a_n_s2_prev=zeros(npos,npos);%a's for spin down at each N at previous time
a_n_s3_prev=zeros(npos,npos);%a's for spin down at each N at previous time
a_n_s4_prev=zeros(npos,npos);%a's for spin down at each N at previous time

a_n_s1_cur=zeros(npos,npos);%a's for spin up at each N at current time
a_n_s2_cur=zeros(npos,npos);%a's for spin down at each N at current time
a_n_s3_cur=zeros(npos,npos);%a's for spin down at each N at current time
a_n_s4_cur=zeros(npos,npos);%a's for spin down at each N at current time

pos_n_counter=0;
pos_m_counter=0;
a=0;
p_n_tmax = zeros(npos,npos);
% p_n = zeros(npos,npos);
x_track=[nmid_rel];
y_track=[nmid_rel];

for iter=1:n_iterations
    display(iter);
    a = 1/2;
%     a = rand;
    a_n_s1_prev(nmid_rel,nmid_rel)=a;
    a_n_s2_prev(nmid_rel,nmid_rel)=a*i;
    a_n_s3_prev(nmid_rel,nmid_rel)=a*i;
    a_n_s4_prev(nmid_rel,nmid_rel)=a;
    
    
    %for all times
    for t=time
        display(t);
        nstart_rel = nmid_rel - t;
        nend_rel = nmid_rel + t;
        
        for pos_n_counter=nstart_rel:nend_rel
%             display(pos_n_counter);
            for pos_m_counter=nstart_rel:nend_rel
%                   display(pos_m_counter);
                %%%% Hadamard Coin
                a_n_s4_cur(pos_n_counter,pos_m_counter) = 0.5 * (a_n_s1_prev( pos_n_counter-1,pos_m_counter-1) + a_n_s2_prev( pos_n_counter-1,pos_m_counter-1) + a_n_s3_prev( pos_n_counter-1,pos_m_counter-1) + a_n_s4_prev( pos_n_counter-1,pos_m_counter-1) );
                a_n_s2_cur(pos_n_counter,pos_m_counter) = 0.5 * (a_n_s1_prev( pos_n_counter-1,pos_m_counter+1) - a_n_s2_prev( pos_n_counter-1,pos_m_counter+1) + a_n_s3_prev( pos_n_counter-1,pos_m_counter+1) - a_n_s4_prev( pos_n_counter-1,pos_m_counter+1) );
                a_n_s3_cur(pos_n_counter,pos_m_counter) = 0.5 * (a_n_s1_prev( pos_n_counter+1,pos_m_counter-1) + a_n_s2_prev( pos_n_counter+1,pos_m_counter-1) - a_n_s3_prev( pos_n_counter+1,pos_m_counter-1) - a_n_s4_prev( pos_n_counter+1,pos_m_counter-1) );
                a_n_s1_cur(pos_n_counter,pos_m_counter) = 0.5 * (a_n_s1_prev( pos_n_counter-1,pos_m_counter-1) - a_n_s2_prev( pos_n_counter+1,pos_m_counter+1) - a_n_s3_prev( pos_n_counter+1,pos_m_counter+1) + a_n_s4_prev( pos_n_counter+1,pos_m_counter+1) );
                     
                %%%Grover coin
%                 a_n_s1_cur(pos_n_counter,pos_m_counter) = 0.5 * (-a_n_s1_prev( pos_n_counter-1,pos_m_counter-1) + a_n_s2_prev( pos_n_counter-1,pos_m_counter-1) + a_n_s3_prev( pos_n_counter-1,pos_m_counter-1) + a_n_s4_prev( pos_n_counter-1,pos_m_counter-1) );
%                 a_n_s2_cur(pos_n_counter,pos_m_counter) = 0.5 * (a_n_s1_prev( pos_n_counter-1,pos_m_counter+1) - a_n_s2_prev( pos_n_counter-1,pos_m_counter+1) + a_n_s3_prev( pos_n_counter-1,pos_m_counter+1) + a_n_s4_prev( pos_n_counter-1,pos_m_counter+1) );
%                 a_n_s3_cur(pos_n_counter,pos_m_counter) = 0.5 * (a_n_s1_prev( pos_n_counter+1,pos_m_counter-1) + a_n_s2_prev( pos_n_counter+1,pos_m_counter-1) - a_n_s3_prev( pos_n_counter+1,pos_m_counter-1) + a_n_s4_prev( pos_n_counter+1,pos_m_counter-1) );
%                 a_n_s4_cur(pos_n_counter,pos_m_counter) = 0.5 * (a_n_s1_prev( pos_n_counter-1,pos_m_counter-1) + a_n_s2_prev( pos_n_counter+1,pos_m_counter+1) + a_n_s3_prev( pos_n_counter+1,pos_m_counter+1) - a_n_s4_prev( pos_n_counter+1,pos_m_counter+1) );
                
            end;
        end;
        
        %record the probability densities of up,down and all particles
        p_n_tmax =  abs(a_n_s1_cur).^2 + abs(a_n_s2_cur).^2  + abs(a_n_s3_cur).^2  + abs(a_n_s4_cur).^2;     
        
%         p_n = p_n + p_n_tmax;
        %set the a arrays of prev time with values of that of the current
        %time and a arrays of current time with zero values
%         if(t <tmax)
            a_n_s1_prev = a_n_s1_cur   ;
            a_n_s2_prev = a_n_s2_cur   ;
            a_n_s3_prev = a_n_s3_cur   ;
            a_n_s4_prev = a_n_s4_cur   ;
            
            a_n_s1_cur=zeros(npos,npos);%a's for spin up at each N at current time
            a_n_s2_cur=zeros(npos,npos);%a's for spin down at each N at current time
            a_n_s3_cur=zeros(npos,npos);%a's for spin down at each N at current time
            a_n_s4_cur=zeros(npos,npos);%a's for spin down at each N at current time
%         end;
        
        if(iter ~= n_iterations) 
            p_n_tmax = zeros(npos,npos);
        end;
        
    end;
   
end;

% p_n = p_n/(n_iterations*tmax);


figure(10)
mesh(pos-nmid_rel,pos-nmid_rel,p_n_tmax);
title('plot of P_n(t) for 2D Quantum hadamard random walk');
xlabel('X axis');
ylabel('Y axis');
grid on;